MODULE inout
  
  USE definitions
  
  IMPLICIT NONE
  
CONTAINS
  
  SUBROUTINE output(i,j)
    
    USE definitions
    USE functions
    
    IMPLICIT NONE
    
    INTEGER, INTENT(in) :: i,j
    INTEGER :: k, x, y, covered_slabs
    REAL :: d, m
    REAL :: avhight, averhight, avoutflow, avvolflux, avtrlength
    IF (mod(j,xsize*ysize)==0 .AND. j /= 0) THEN
       d=0
       DO k=0,5000
          d = d + path(k)*k
          path(k)=0
       ENDDO

       CALL masscount(avhight)
       averhight = erbuf/real(xsize*ysize)
       avoutflow = xsinkbuf/real(ysize)
       avtrlength = d/(REAL(pn)+1)
       avvolflux = averhight*avtrlength
       
       erbuf = 0
       xsinkbuf = 0

       IF(avhightbuf>avhight) THEN
          qtrig = 1
          IF(tconv == 0) tconv = i
       ENDIF
       WRITE(*,*) qtrig, avhightbuf, avhight, avhight-avhightbuf
       avhightbuf = avhight
       
       WRITE(*,*) 'Progress:', i, 'of', steps, 'steps', char(9)
       WRITE(*,*) d/(REAL(pn)+1), 'Mean path length'
       WRITE(*,*) avhight, 'Average hight per Cell'
       WRITE(*,*) avoutflow, 'Mass input/output per cell on boundary'
       WRITE(*,*) avvolflux, 'average volume flux on grid'
       WRITE(out4, '(I5, 7E16.8)') i, avtrlength, avhight, averhight,&
            & avoutflow, avvolflux, dwflux, daflux
       WRITE(*,*) 'qtrig=', qtrig, 'tconv=', tconv
       empty = 0
       pn = 0
       
       IF(slices == 1) CALL slice_diagnostics(i)
       covered_slabs = 0
       DO y=1,ysize
          DO x=1,xsize
             IF(sandpit(x,y)%slabnumber > max_hight(x))max_hight(x) =&
                  & sandpit(x,y)%slabnumber
             IF(sandpit(x,y)%slabnumber > local_max_hight(x,y))&
                  & local_max_hight(x,y) = sandpit(x,y)%slabnumber
             IF(sandpit(x,y)%slabnumber > 0) covered_slabs =&
                  & covered_slabs + 1
          ENDDO
       ENDDO
       
    ENDIF
    
    IF (mod(j,xsize*ysize)==0 .AND. mod(i,frames) == 0) THEN
       CALL write_slabnumber(i/frames) !Write current system state to
       ! file
       CALL SYSTEM('sh mvframe.sh')
       IF(full == 1) CALL full_diagnostics(i/frames)
       IF(full == 1) CALL SYSTEM('sh mvdiag.sh')
    ENDIF

       IF(full == 1 .AND. mod(j,xsize*ysize)==0 .AND. j /= 0) THEN

          local_transport_length = 0
          local_eroded_hight = 0
          local_deposited_hight = 0
          local_aeolian_flux = 0
          local_avalanche_flux = 0
          local_sediment_balance = 0
       ENDIF

  END SUBROUTINE output
  !------------------------------------------------------------

  !------------------------------------------------------------
  !Run slice diags
  !------------------------------------------------------------
  SUBROUTINE slice_diagnostics(step)

    USE definitions
    USE functions
    USE transport
    USE erosion

    IMPLICIT NONE

    INTEGER :: x, y
    INTEGER, DIMENSION(2) :: pos
    REAL :: L, buff, tmp
    INTEGER, INTENT(in) :: step
    REAL, DIMENSION(:), ALLOCATABLE :: slice
    REAL :: avtrlength, averhight, avhight, avflux, covered_cells

   allocate(slice(1:xsize))

    DO y=1,ysize
       
       slice = sandpit(:,y)%slabnumber

       avtrlength = 0
       averhight = 0 
       avhight = 0
       avflux = 0
       covered_cells = 0

       DO x=1,xsize
          tmp = sandpit(x,y)%slabnumber
          pos(1) = x
          pos(2) = y

          CALL calcl(pos, L)
          avtrlength = avtrlength + L*(1+par2*ATAN(gradient_h(pos)))
          averhight = averhight + erode(pos, 0)
          avhight = avhight + sandpit(x,y)%slabnumber
          pos(1) = pos(1) + INT(2.0*L0/3.0)
          IF(pos(1) < xsize) avflux = avflux + L*(1+par2*ATAN(gradient_h(pos)))*erode(pos, 0)
          IF(sandpit(x,y)%slabnumber > 0) covered_cells = covered_cells + 1
          IF(tmp /= sandpit(x,y)%slabnumber) write(*,*) 'Danger, Danger!'
       ENDDO

       sliced_transport_length(step,y) = avtrlength/xsize
       sliced_erosion_hight(step,y) = averhight/xsize 
       sliced_average_hight(step,y) = avhight/xsize
       sliced_average_flux(step,y) = avflux/xsize
       sliced_covered_cells(step,y) = covered_cells/xsize

    ENDDO

    deallocate(slice)

  END SUBROUTINE slice_diagnostics
  !------------------------------------------------------------

  !------------------------------------------------------------
  !routine for full diagnostics including local sediment balance, 
  !local transport flux, eroded hight and eroded volume. output 
  !is written to different files for each quantity and each 
  !designated time step
  !------------------------------------------------------------
  SUBROUTINE full_diagnostics(step)

    USE definitions
    USE transport
    USE erosion
    USE deposition
    USE functions
    
    IMPLICIT NONE
    
    INTEGER :: irec, step, x, y, counter
    CHARACTER(len=31) :: filename
    CHARACTER(len=5) :: ci
!Collect stuff for diagnostics and format it for output

    local_sediment_balance = local_deposited_hight - local_eroded_hight

!Write arrays to files according to name and step number    
    WRITE(ci, '(I5.5)') step
    ci = TRIM(ADJUSTL(ci))
    filename = 'local_transport_length' // ci // '.dat'

    INQUIRE(iolength = irec) local_transport_length
    OPEN(unit=out6, file=filename, form = 'UNFORMATTED', access='direct', recl=irec )
    WRITE(out6, rec = 1) local_transport_length
    CLOSE(out6, status = 'keep')

    WRITE(ci, '(I5.5)') step
    ci = TRIM(ADJUSTL(ci))
    filename = 'local_eroded_hight' // ci // '.dat'

    INQUIRE(iolength = irec) local_eroded_hight
    OPEN(unit=out7, file=filename, form = 'UNFORMATTED', access='direct', recl=irec )
    WRITE(out7, rec = 1) local_eroded_hight
    CLOSE(out7, status = 'keep')

    WRITE(ci, '(I5.5)') step
    ci = TRIM(ADJUSTL(ci))
    filename = 'local_deposited_hight' // ci // '.dat'

    INQUIRE(iolength = irec) local_deposited_hight
    OPEN(unit=out7, file=filename, form = 'UNFORMATTED', access='direct', recl=irec )
    WRITE(out7, rec = 1) local_deposited_hight
    CLOSE(out7, status = 'keep')

    WRITE(ci, '(I5.5)') step
    ci = TRIM(ADJUSTL(ci))
    filename = 'local_aeolian_flux' // ci // '.dat'

    INQUIRE(iolength = irec) local_aeolian_flux
    OPEN(unit=out8, file=filename, form = 'UNFORMATTED', access='direct', recl=irec )
    WRITE(out8, rec = 1) local_aeolian_flux
    CLOSE(out8, status = 'keep')

    WRITE(ci, '(I5.5)') step
    ci = TRIM(ADJUSTL(ci))
    filename = 'local_avalanche_flux' // ci // '.dat'

    INQUIRE(iolength = irec) local_avalanche_flux
    OPEN(unit=out8, file=filename, form = 'UNFORMATTED', access='direct', recl=irec )
    WRITE(out8, rec = 1) local_avalanche_flux
    CLOSE(out8, status = 'keep')

    WRITE(ci, '(I5.5)') step
    ci = TRIM(ADJUSTL(ci))
    filename = 'local_sediment_balance' // ci // '.dat'

    INQUIRE(iolength = irec) local_sediment_balance
    OPEN(unit=out9, file=filename, form = 'UNFORMATTED', access='direct', recl=irec )
    WRITE(out9, rec = 1) local_sediment_balance
    CLOSE(out9, status = 'keep')

    mintl=local_transport_length(1,1)
    maxtl=local_transport_length(1,1)
    mineh=local_eroded_hight(1,1)
    maxeh=local_eroded_hight(1,1)
    mindh=local_deposited_hight(1,1)
    maxdh=local_deposited_hight(1,1)
    minaf=local_aeolian_flux(1,1)
    maxaf=local_aeolian_flux(1,1)
    minvf=local_avalanche_flux(1,1)
    maxvf=local_avalanche_flux(1,1)
    minsb=local_sediment_balance(1,1)
    maxsb=local_sediment_balance(13,1)

    dwflux = 0
    daflux = 0
    counter = 0
    DO y=1,ysize
       DO x=1,xsize
          IF(sandpit(x,y)%slabnumber>2.0*E0) THEN
             counter = counter+1
             dwflux = dwflux + local_avalanche_flux(x,y)
             daflux = daflux + local_aeolian_flux(x,y)
          ENDIF
          IF(local_transport_length(x,y)<mintl) mintl=local_transport_length(x,y)
          IF(local_transport_length(x,y)>maxtl) maxtl=local_transport_length(x,y)
          IF(local_eroded_hight(x,y)<mineh) mineh=local_eroded_hight(x,y)
          IF(local_eroded_hight(x,y)>maxeh) maxeh=local_eroded_hight(x,y)
          IF(local_deposited_hight(x,y)<mindh) mindh =local_deposited_hight(x,y)
          IF(local_deposited_hight(x,y)>maxdh) maxdh=local_deposited_hight(x,y)
          IF(local_aeolian_flux(x,y)<minaf) minaf=local_aeolian_flux(x,y)
          IF(local_aeolian_flux(x,y)>maxaf) maxaf=local_aeolian_flux(x,y)
          IF(local_avalanche_flux(x,y)<minvf) minvf=local_avalanche_flux(x,y)
          IF(local_avalanche_flux(x,y)>maxvf) maxvf=local_avalanche_flux(x,y)
          IF(local_sediment_balance(x,y)<minsb) minsb=local_sediment_balance(x,y)
          IF(local_sediment_balance(x,y)>maxsb) maxsb=local_sediment_balance(x,y)
       ENDDO
    ENDDO

    dwflux = dwflux/REAL(counter)
    daflux = daflux/REAL(counter)

    local_transport_length = 0
    local_eroded_hight = 0
    local_deposited_hight = 0
    local_aeolian_flux = 0
    local_avalanche_flux = 0
    local_sediment_balance = 0

  END SUBROUTINE full_diagnostics
  !------------------------------------------------------------
  
  !------------------------------------------------------------
  !Function to load simulation input parameters from ast.in
  !------------------------------------------------------------
  SUBROUTINE load_parameters
    
    USE definitions
    
    INTEGER, PARAMETER :: in=10
    
    OPEN(unit=in, file="ast.in")
    READ(in, PARAMETER)
    CLOSE(in)
    
  END SUBROUTINE load_parameters
  !------------------------------------------------------------
  
  
  
  !------------------------------------------------------------
  !Function to write slab hight profile to hight.dat
  !------------------------------------------------------------
  SUBROUTINE write_slabnumber(num)
    
    USE definitions
    USE functions
    USE deposition
    USE erosion
    
    IMPLICIT NONE
    
    INTEGER 			:: i, j, num, ls1, ls2, k, irec
    INTEGER, DIMENSION(2) :: pos
    CHARACTER(len=5) :: ci
    CHARACTER(len=15) :: filename, str1, str2=''
    REAL, DIMENSION(xsize, ysize) :: tmp
    !make string for output file name
    
    WRITE(ci, '(I5.5)') num
    ci = TRIM(ADJUSTL(ci))
    filename = 'hight' // ci // '.dat'
    
    !Write data to filename.dat

       DO j=1, ysize	
          DO i=1, xsize
             tmp(i,j) =  REAL(Sandpit(i,j)%Slabnumber)
          END DO
       END DO

       INQUIRE(iolength = irec) tmp
       OPEN(unit=out1, file=filename, form = 'UNFORMATTED', access='direct', recl=irec )
       WRITE(out1, rec = 1) tmp
    CLOSE(out1, status = 'keep')

  END SUBROUTINE write_slabnumber
  !------------------------------------------------------------
  !Routine to write final state of the System to file for continuing simulation
  !------------------------------------------------------------
  SUBROUTINE write_final_state
    
    USE definitions
    USE functions
    USE deposition
    USE erosion
    
    IMPLICIT NONE
    
    INTEGER :: i, j, num, ls1, ls2
    INTEGER, DIMENSION(2) :: pos
    !Write data to final.dat
    
    OPEN(unit=out2, file='final.dat')
    
    WRITE(out2, '(2I4)') xsize, ysize
    
    DO i=1, xsize
       DO j=1, ysize
          
          pos(1)=i
          pos(2)=j
          
          WRITE(out2, '(1I5,A, 1I5, A, E16.8)') &
               i, char(9), j, char(9), Sandpit(i,j)%Slabnumber
          
       END DO
    END DO
    
    CLOSE(unit = out2)

  END SUBROUTINE write_final_state
  !------------------------------------------------------------
  
  
  !------------------------------------------------------------
  !Routine to write simulation Parameters to file for readout via bash/gnuplot
  !------------------------------------------------------------
  SUBROUTINE write_parameters
    
    USE definitions
    
    IMPLICIT NONE
    
    
    INTEGER :: i, j
    REAL    :: minhight=0, maxhight=0, maxshadow=0, minshadow=0
    
    minhight = REAL(initialhight)
    
    DO j=1,ysize
       DO i=1,xsize
          IF(sandpit(i,j)%slabnumber>maxhight) maxhight = sandpit(i,j)%slabnumber
          IF(sandpit(i,j)%slabnumber<minhight) minhight = sandpit(i,j)%slabnumber
       ENDDO
    ENDDO

    OPEN(unit=out4, file = 'parameters.dat')
  
    WRITE(out4, '(A, I4, A, A5)') '#', xsize, char(9), 'xsize'
    WRITE(out4, '(A, I4, A, A5)') '#', ysize, char(9), 'ysize'
    WRITE(out4, '(A, I4, A, A5)') '#', steps+1, char(9), 'steps'
    WRITE(out4, '(A, E16.8, A, A10)') '#', sourceflow, char(9), 'sourceflow'
    WRITE(out4, '(A, E16.8, A, A2)') '#', L0, char(9), 'L0'
    WRITE(out4, '(A, E16.8, A, A2)') '#', E0, char(9), 'E0'
    WRITE(out4, '(A, E16.8, A, A8)') '#', minhight, char(9), 'minhight'
    WRITE(out4, '(A, E16.8, A, A8)') '#', maxhight, char(9), 'maxhight'
    WRITE(out4, '(A, E16.8, A, A8)') '#', slabsize, char(9), 'slabsize'
    WRITE(out4, '(A, I4, A, A5)') '#', tconv, char(9), 'tconv' 
    
    CLOSE(unit = out4)
    
  END SUBROUTINE write_parameters

  SUBROUTINE write_envelope

    USE definitions

    IMPLICIT NONE

    INTEGER :: y, irec

    OPEN(unit = out20, file = '1d_envelope.dat')

    DO y=1,ysize
       WRITE(out20, '(I4, E16.8)') y, max_hight(y)
    ENDDO

    CLOSE(unit = out20)

    INQUIRE(iolength = irec) local_max_hight
    OPEN(unit = out20, file = '2d_envelope.dat', form = 'UNFORMATTED', access='direct', recl=irec)

    WRITE(unit = out20, rec = 1) local_max_hight

    CLOSE(unit = out20)

  END SUBROUTINE write_envelope
  !------------------------------------------------------------
  
  SUBROUTINE fft
    
    USE definitions

    IMPLICIT NONE
    
    integer ( kind = 4 ) :: l
    integer ( kind = 4 ) :: m
    integer ( kind = 4 ) :: irec
    complex ( kind = 8 ), allocatable :: c(:,:)
    real, allocatable :: cout(:,:)
    integer ( kind = 4 ) ier
    integer ( kind = 4 ) ldim
    integer ( kind = 4 ) lensav
    integer ( kind = 4 ) lenwrk
    real ( kind = 8 ), allocatable, dimension ( : ) :: work
    real ( kind = 8 ), allocatable, dimension ( : ) :: wsave

    WRITE(*,*) 'Calculating fourier diagnostics of the hight profile'
    
    l = xsize
    m = ysize
    
    allocate (    c(xsize,ysize) )
    allocate ( cout(xsize,ysize) )
    
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TEST17'
    write ( *, '(a)' ) '  For complex double precision fast Fourier transforms, 2D,'
    write ( *, '(a)' ) '  ZFFT2I initializes the transform,'
    write ( *, '(a)' ) '  ZCFFT2F does a forward transform;'
    write ( *, '(a)' ) '  ZFFT2B does a backward transform.'
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  The data is stored in an L by M array, with'
    write ( *, '(a,i8)' ) '  L = ', l
    write ( *, '(a,i8)' ) '  M = ', m
    
    !
    !  Set work arrays.
    !
    lensav = 2 * ( l + m ) + int ( log ( real ( l, kind = 8 ) ) ) &
         + int ( log ( real ( m, kind = 8 ) ) ) + 8
    lenwrk = 2 * l * m
    
    write ( *, '(a,i8)' ) '  LENSAV = ', lensav
    write ( *, '(a,i8)' ) '  LENWRK = ', lenwrk
    
    allocate ( wsave(1:lensav) )
    allocate ( work(1:lenwrk) )
    
    call zfft2i ( l, m, wsave, lensav, ier )
    
    !
    !  Set the data values.
    !
    
    c = CMPLX(Sandpit(:,:)%Slabnumber)
    
    !
    !  Compute the FFT coefficients.
    !
    
    ldim = l
    
    call zfft2f ( ldim, l, m, c, wsave, lensav, work, lenwrk, ier )
    
    cout = REAL(c)

    INQUIRE(iolength = irec) cout
    OPEN(unit=out5, file ='fouriercomps.dat', form = 'UNFORMATTED', access='direct', recl=irec)
    WRITE(unit = out5, rec = 1) cout
    CLOSE(unit = out5)
    
    deallocate ( work )
    deallocate ( wsave )
    
    return
    
  END SUBROUTINE fft
 END MODULE inout
